

class TransdimGibbsst
{
	public:
		TransdimGibbsst(mat v,mat X,double *Y, double ds)
		{
			_ds=ds;
			_v=v;
			rt=new rtnorm();
			_X=X;
			_p=X.n_cols-1;
			_n=X.n_rows;
			_Y=Y;
			U=new Distribution::Uniform(1,0,1);
			boost::random::uniform_01<> U;
			runif=new RandomG::Random<boost::random::uniform_01<> >(U);
			_prec=10;
			_vp=diagmat(v);
		}
		mat Pi_b(mat gamma, mat z,mat b)
		{
			mat h=sum(gamma,0);
			int nb=h(0,0)+1;
			mat gammat(_p+1,1);
			mat X2(_n,nb);
//			cout << "test"<< p+1;
			int kk=1;
			for(int i=0;i<(_p+1);i++)
			{
					
				if(i!=0){
					gammat(i,0)=gamma(i-1,0);
					if(gamma(i-1,0)!=0)
					{
						X2(span::all,kk)=_X(span::all,i);
						kk++;
					}

				}else{
					gammat(i,0)=1;
					X2(span::all,0)=_X(span::all,0);
				}
			}
			mat v=Adjust(_vp,gammat);
			mat V=inv(inv(v)+X2.t()*X2);
			mat m=V*X2.t()*z;	
			/*cout << V;
			cout << X2;
			cout << nb;*/
			G=new Distribution::Gaussian(nb,m,V);
			mat res1=(*G).r(1);
			mat beta=add_mat(res1,gammat,_p+1);
			//cout << z.t();
			//cout << beta.t();	
			delete G;	
			return beta;			
		}
		mat Q_gamma(mat gammat)
		{
			Distribution::Uniform U1(1,0,_p-1);
			int u=as_scalar(U1.r(1));
			gammat(u,0)=1-gammat(u,0);//as_scalar((*U).r(1));	
			return gammat;
		}
		void rtau(mat b)
		{
			for(int i=0;i<_p;i++)
			{
				double s=(pow(_prec,2)+pow(b(0,i),2))/2;
				boost::random::gamma_distribution<> Ni(1,1/s);
				RandomG::Random<boost::random::gamma_distribution<> > rgamma(Ni);
				_vp(i,i)=1/(rgamma)();
			}
		}
		vector<mat> prob(mat gamma,mat gammat1,mat z,mat b,mat bt)
		{
			mat h=sum(gamma,0);
			int nb=h(0,0)+1;
			vector<mat> vect;
			mat gammat(_p+1,1);
			mat gammat1t(_p+1,1);
			mat X2(_n,nb);
//			cout << "test"<< p+1;
			int kk=1;
			for(int i=0;i<(_p+1);i++)
			{
					
				if(i!=0){
					gammat(i,0)=gamma(i-1,0);
					gammat1t(i,0)=gammat1(i-1,0);
					
					if(gamma(i-1,0)!=0)
					{
						X2(span::all,kk)=_X(span::all,i);
						kk++;
					}

				}else{
					X2(span::all,0)=_X(span::all,0);
					gammat(i,0)=1;
					gammat1t(i,0)=1;
				}
			}
			
			mat v=Adjust(_vp,gammat);
			mat V=inv(inv(v)+X2.t()*X2);
			mat B=V*X2.t()*z;
			mat Vg=inv(inv(_v)+_X2.t()*_X2);
			mat Bg=Vg*_X2.t()*z;
		
			double num=as_scalar(0.5*log(det(V))+0.5*log(det(_v))+(0.5*B.t()*inv(V)*B)+Prior(gammat));
			double denom=as_scalar(0.5*log(det(Vg))+0.5*log(det(v))+(0.5*Bg.t()*inv(Vg)*Bg)+Prior(gammat1t));
			double res=num-denom;
			//cout << "res: "<< res;
			double u=log((*runif)());
			if(u<res)
			{
				_X2=X2;
				_v=v;
				vect.push_back(b);
				vect.push_back(gamma);

			}else{
				vect.push_back(bt);
				vect.push_back(gammat1);
			}
			return vect;
		}

		vector<mat> P_gamma_b(mat gammat,mat z, mat bt)
		{
			mat gamma=Q_gamma(gammat);
			mat b=Pi_b(gamma,z,bt);
			vector<mat> res=prob(gamma,gammat,z,b,bt);
			return res;
		}
		mat rz(mat b)
		{
			mat z(_n,1);
			for(int i=0;i<_n;i++)
			{
				double foo=dot(_X.row(i),b);
				if(_Y[i]==0)
				{
					z(i,0)=(*rt)(foo,-pINF,0,1);
				}else{
					z(i,0)=(*rt)(foo,0,pINF,1);
				}
			//	cout << z(i,0) << " y " << _Y[i] << "\n"; 
			}
			return z;
		}
		mat Sample(int M)
		{
			mat Gamma(M,_p);
			Gamma.fill(0);
			post_init(Gamma.row(0).t());	
			mat z=_z;
			mat b(_p+1,1);
			b.fill(0);
			b(0,0)=_Bg(0,0);
			for(int i=1;i<M;i++)
			{
				vector<mat> vect=P_gamma_b((Gamma.row(i-1)).t(),z,b);
				b=vect[0];
				Gamma.row(i)=vect[1].t();
				cout << Gamma.row(i);
				z=rz(b);
				rtau(vect[0].t());
				//cout << b.t();
				cout << " " << i;
			}
			return Gamma;

		}

		void post_init(mat gamma)
		{

			int nkp=as_scalar(sum(gamma,0));
			int nb=nkp+1;
			vector<mat> vect;
			vector<mat> res;
			mat gammati(_p+1,1);
			mat X2(_n,nb);
//			cout << "test"<< p+1;
			int kk=1;
			
			for(int i=0;i<(_p+1);i++)
			{
					
				if(i!=0){
					gammati(i,0)=gamma(i-1,0);
					if(gamma(i-1,0)!=0)
					{
						X2(span::all,kk)=_X(span::all,i);
						kk++;
					}

				}else{
					X2(span::all,0)=_X(span::all,0);
					gammati(i,0)=1;
				}
			}


			mat Xt=_X*diagmat(gammati);
			mat v=Adjust(_vp,gammati);
			_X2=X2;
			_Vg=inv(inv(v)+X2.t()*X2);
			mat Y=vect2mat<double>(_Y,_n,1);
			_v=v;
			_Bg=_Vg*X2.t()*Y;
			_z=X2*_Bg;	

		}

		double Prior(mat theta){
		    
			double t=0;
			double m=0.5;
			double ds=logfact(_ds);
			double d=logfact(_p);
			double i=sum(sum(theta));
			if(i<=_ds)
			{
				double t=logfact(i)+(logfact(_ds-i))-ds+i*log(m)+(_ds-i)*log(1-m)+logfact(i)+(logfact(_p-i))-d;
				return t;
 			}else{
				return -numeric_limits<double>::infinity();

			}		
			
	//		return -0.693;
		}



	private:
		rtnorm *rt;
		Distribution::Gaussian *G;
		mat _v;
		Distribution::Uniform *U;
		double *_Y;
		mat _X;
		int _p;
		int _n;
		int _ds;
		double _prec;
		RandomG::Random<boost::random::uniform_01<> > *runif;
		mat _Bg;
		mat _X2;
		mat _Vg;
		mat _z;
		mat _vp;

};
